Post MCMCglmm analyses

In this script I’ll be looking at how to analyse the outputs from the MCMCglmm models to get the something about the elaboration/exploration stories.

I’m going to look at it using two approaches: * the blob-wise one where I’ll look at the differences between whole blobs (ellipses) in terms of elaboration/exploration (Robinson & Beckerman style); * and the tip-wise one where I’ll look at the elaboration/exploration score of the tips on the tree (cf. their blobs) relative to the whole phylogeny and to their respective blobs (clades).

Data

I’m going to focus on the data from Gavin and especially the models 5 and 6 that are:

• model 6: model_phylo1_clade3 with one phylogenetic random effect (across the whole tree) and three clade residual effects.
• model 7: model_phylo3_clade3 with three phylogenetic random effect (one for each clade) and three clade residual effects.

## Loading the correct models
load("../Data/Processed/model_list.rda")
model_phylo1_clade3 <- model_list[[4]]
model_phylo3_clade3 <- model_list[[6]]

And here’s what the trait space looks like:

## Loading the data
load("../Data/Processed/morphdat.rda")
## Plotting it
colour_vector <- c("orange", "blue", "darkgreen")
plot(morphdat[, c(1,2)], pch = 19,
     col = colour_vector[morphdat$clade],
     xlab = "PC1 (90.5%)",
     ylab = "PC1 (6.86%)")
legend("topleft", legend = levels(morphdat$clade), col = colour_vector, pch = 19)

Note that the PC% are from the whole PC in Gavin’s example.

Blob-wise approach

Here we’re basically comparing multidimensional ellipses (here 3D ones but I think we should push it to more dimensions!).

Phylo + clade model

First we can visualise some of these ellipses (100, randomly drawn from the posterior):

source("../Functions/plot.ellipses.R")
source("../Functions/get.covar.R")
plot(morphdat[, c(1,2)], pch = 19, cex = 0.5,
     col = colour_vector[morphdat$clade],
     xlab = "PC1 (90.5%)",
     ylab = "PC1 (6.86%)")
plot.ellipses(model_phylo1_clade3, n = 100,
              centre = "none", add = TRUE, col = c("grey",colour_vector))

Because this is really messy, we can centre these ellipses on the groups ellipses average centres:

plot(morphdat[, c(1,2)], pch = 19, cex = 0.5,
     col = colour_vector[morphdat$clade],
     xlab = "PC1 (90.5%)",
     ylab = "PC1 (6.86%)")
plot.ellipses(model_phylo1_clade3, n = 100,
              centre = "level", add = TRUE, col = c("grey",colour_vector))

Phylo clade + clade model

plot(morphdat[, c(1,2)], pch = 19, cex = 0.5,
     col = colour_vector[morphdat$clade],
     xlab = "PC1 (90.5%)",
     ylab = "PC1 (6.86%)")
plot.ellipses(model_phylo3_clade3, n = 100,
              centre = "none", add = TRUE, col = c("black", colour_vector, colour_vector))

This looks much more messier (and also probably not scaled correctly?)…

plot(morphdat[, c(1,2)], pch = 19, cex = 0.5,
     col = colour_vector[morphdat$clade],
     xlab = "PC1 (90.5%)",
     ylab = "PC1 (6.86%)")
plot.ellipses(model_phylo3_clade3, n = 100,
              centre = "level", add = TRUE, col = c("black", colour_vector, colour_vector))

Tip-wise approach